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We obtain multi-soliton solutions of the time-dependent Bogoliubov-de Gennes equations or, equivalently, 
Gorkov equations that describe the dynamics of a fermionic condensate in the dissipationless regime. There 



f^ ' are two kinds of sohtons - normal and anomalous. At large times, normal multi-solitons asymptote to 

unstable stationary states of the BCS Hamiltonian with zero order parameter (normal states), while the 
r^ I anomalous ones tend to eigenstates characterized by a nonzero anomalous average. Under certain circum- 

fj ' stances, multi-soliton solutions break up into sums of single solitons. In the linear analysis near the stationary 

5-H ' states, solitons correspond to unstable modes. Generally, they are nonlinear extensions of these modes, so 

^ , that a stationary state with k unstable modes gives rise to a fc-soliton solution. We relate parameters of 



the multi-solitons to those of the asymptotic stationary state, which determines the conditions necessary for 
exciting solitons. We further argue that the dynamics in many physical situations is multi-soliton. 
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INTRODUCTION AND SUMMARY OF THE RESULTS 



Recent years have witnessed renewed experimental and theoretical interest in far from equi- 
librium phenomena in strongly interacting many-body systems at low temperatures. Examples 



include non-stationary Kondo and other impurity models 



W 



14 



£|, quenched Luttinger liquids lCll4l3l|. 



22( 1 etc. On the theory side, there 



electron spin dynamics induced by hyperfine interactions 
is a considerable effort to develop new approaches to nonequilibrium many-body physics. This 
presents a significant challenge as conventional techniques are often inadequate for the description 
of these phenomena. In particular, there have been major advances in the theory of dynamical 
fermionic pairing in the collisionless regime 23143 7l| . This problem is long known to be accurately 



described by the time-dependent Bogoliubov-de Gennes equations, which in this case are a set 

3, |30(]. Nevertheless, it was not until 
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25, m 



29 



301, 



331]. The exact 



of coupled nonlinear integro-differential equations 

recently that this nonlinear system was realized to be exactly solvable [! 
solution proved to be a unique approach to the problem of dynamical pairing and has been exten- 
sively exploited to obtain analytical information about its key physical properties. For example, 
a nonequilibrium "phase diagram" of a homogeneous Bardeen-Cooper-Schriffer (BCS) superfluid 
with a number of novel pha ses , as well as their responses to existing experimental probes were 
predicted analytically 



34 



I35|,l3a,l37(. 

However, while much attention was focused on the asymptotic states of the condensate at large 
times, the transient dynamics has not been fully explored. Most importantly, nonlinear integrable 
systems are known to exhibit a remarkable class of multi-soliton solutions that play a central role 
in understanding and predicting their properties. Physical solutions can often be represented as a 
superposition of solitons making a quantitative analysis possible. For instance, one can show that 
the dynamics giving rise to nonequilibrium "phases" mentioned above is multi-soliton, see below. 
The existence and properties of solutions of this type are also of a general interests from the point 
of view of nonlinear physics due to the nonlocal nature of the BCS problem distinguishing it from 
familiar integrable systems, such as nonlinear Shrodinger, Korteweg-de Vries, sine-Gordon etc. 

In this paper we construct multi-soliton solutions to the dynamical fermionic pairing problem, 
see Figs. [D - El We establish a direct correspondence between solitons and the stationary states 
of the mean-field BCS Hamiltonian, such that each soliton solution asymptotes to an eigenstate 
at times t — > itcxD. This also identifies conditions necessary for exciting solitons. There are two 
distinct types of solitons - normal and anomalous. Normal solitons asymptote to stationary states 
that are simultaneous eigenstates of a Fermi gas and the mean-field BCS Hamiltonian and are 



characterized by a zero anomalous average. For anomalous solitons the asymptotic value of this 
average is finite, A(t -^ iLcxd) ^ 0. As the separation between solitons is increased, the multi- 
soliton solution breaks up into a simple sum of single solitons (see e.g. Figs. [2] and [3])- one of the 
defining properties of solitons. 

In the rest of this section, we briefly formulate the problem and then summarize the main 
results. The collisionless dynamics of a fermionic superfluid can be described by the Bogoliubov-de 
Gennes equations [27|, |30(] 




^m 




... (1) 

A*(t) -Er. 



where A(t) = g'}2m^rnym is the anomalous average, e^ are the single fermion energies relative 



to the Fermi level ep, and g is the coupling constant. These non 
integrable for any number of Bogoliubov amplitudes {Um, Kn)l29 



inear equations are known to be 



3d, 



33] • In the continuum limit 



the summation in the expression for A(t) is replaced by integration and Eqs. ([T]) become integrable 
nonlinear integro-differential equations. Each solution of Eq. ([T]) yields a many-body wave function 

i^w) = WPmit) + v:;,{t)ci^ci^m, (2) 

m 

where the product is taken only over unblocked levels - levels that are either unoccupied or doubly 
occupied. 

As we demonstrate in subsequent sections, solitons tend to unstable stationary states of the 
mean-field BCS Hamiltonian at t ^ ztcxo. The mean-field Hamiltonian has two types of eigenstates 
- normal and anomalous. Anomalous stationary states have a nonzero constant value A^ of the 



anomalous average and are solutions of Eq. ([T]) of the form|38l. l39l | {Um,Vm) = {Um,Vm)e~^^"^*, 
where {U^, V^) and Em are the eigenvectors and eigenvalues of the 2x2 matrix on the right hand 
side of Eq. ([1]). There are two states Em = ±y^e^ + A^ for each e^- The BCS ground state has 
Em < for all m. A state where Er > while Em^r < describes a single excited pair [3 a. 



and has energy 2y^e'^ + A^ above the ground state. We note also that an excited pair introduces 
a discontinuity in the average fermion occupation number n{em) = 1 — ^m/Em, since Em changes 
sign at m = r. Normal eigenstates have A = and amplitudes {Um, Vm) equal to either (0, e*'^'"*) or 
(e""™*, 0). Both types of eigenstates are exact stationary states of the mean-field BCS Hamiltonian 

where A = 5 X]fc(cA:iCfc|), and &■ and Cj^ are the creation and annihilation operators for the two 
fermion species. Normal states are also eigenstates of the Fermi gas - the first term in Eq. ([3]). For 



example, the Fermi ground state is a normal eigenstate with Um = for Cm < and Vm = for 
Cm > 0, i.e. all single particle states below the Fermi level occupied and states above it empty. 
A linear analysis of Eq. ([1]) around stationary states shows that some of them are unstable 



271, 



41|. These states give rise to solitons, as is typical in integrable nonlinear dynamics. For example, 
a simple pendulum displays a soliton solution when started in its unstable equilibrium with zero 
velocity. In the phase space the soliton is the separatrix connecting the unstable equilibrium 
to itself. The same is true for example for the single soliton solution of the Korteweg-de Vries 
equation 



42l|. Similarly in the present case the unstable modes start to grow exponentially and 
become solitons due to nonlinear effects. In a certain sense, solitons can be viewed as nonlinear 
extensions of the unstable modes. 

Now let us summarize soliton solutions of the Bogoliubov-de Gennes equations ([1]). Detailed 
derivation and discussion of these as well as some other solutions can be found in subsequent 
sections. Here we present only the results for the amplitude of the order parameter |A(t)|. 

A. Normal solitons 

First, we present normal multi-solitons derived in Sec. IIVI The 1-soliton solution of a normal 
type has been previously obtained in Ref 



271 . The amplitude of the order parameter has the 



following form: 



i^«i = JShfey- (^' 



where a is a real parameter. At t ^ ±00 the corresponding wave function ([2]) asymptotes to the 
Fermi ground state. The latter is unstable in the presence of the pairing interaction. Indeed, a 
linear analysis of equations of motion ([T]) around the Fermi ground state shows a single unstable 



normal mode, which grows as |A(t)| oc e '^^ 27|,|4l[ as can be seen from Eq. ^ in the t — > —00 limit 



The plot of Eq. (j3]) is a single peak centered at t = —a/27, see Fig. [TJ Its height (the amplitude 
of the soliton) and width are controlled by the parameter 2j. In the present case 2j = Aq, where 
Aq is the ground state BCS gap. This parameter can be interpreted as an imaginary "frequency" 
of the unstable normal mode, while Eq. (j3|) as an extension of this mode to the nonlinear regime. 
Below we will see that the number of unstable modes and consequently the number of solitons 
corresponding to a given stationary state is related to the number of discontinuities in the average 
fermion occupation number n(e) in this state. Specifically, 2k — 1 discontinuities (jumps) in n(e) 
lead to up to k coupled normal solitons. The Fermi ground state has a single discontinuity at the 
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FIG. 1: (color online) k = 1 normal-soliton solution of the Bogoliubov-de Gennes equations ([T]) 27[, where 
Aq is the ground state BCS gap. At t — > ±cx) the solution asymptotes to the Fermi ground state. The inset 
shows the average fermion occupation number n(e) in this state. A single discontinuity (2k — 1 = 1) in n(e) 
at the Fermi energy gives rise to a single soliton, see the text. 



Fermi level, giving rise to the 1-soliton solution. 

The 2-normal-solitons have considerably richer structure. Let us give two examples. Both are 
described by 

h{t) 



\A{t)\=A 



h{t)h{t)-h?{t) 

with different choices for the amplitude A and function h{t). One choice is A = A\'^2 ~ li\ ^-^d 

^j^^cosh(27it + ai) ^ ^^^ cosh(272t + 0:2 



h{t) 



271 



+ e 



272 



(5) 



(6) 



Examples of |A(t)| for this case are plotted in Fig. [2^. The second option is yl = 16fi\/fi'^ + 7^ 
and 



hit) 



g-2*/.t+i0i cosh(27t + oi - if3) ^ ^2if,t+i4>2 cosh{2jt + 02 + ifi) 



27 



27 



(7) 



where /3 is the phase of /x + 17, i.e. fi + ij = \J ^^ + 7^6*^. Fig. [2)o shows graphs of | A(t)| obtained 
using Eq. ([7]). In both cases the wave function asymptotes to a normal eigenstate at t — > itcxD. This 
eigenstate is obtained from the Fermi ground state by moving all fermions in the energy interval 
—a < em < below the Fermi level to the symmetric interval < e^ < a above it. Eqs. Q and ([7]) 
correspond to different values of a. Note that this normal eigenstate has 2A; — 1 = 3 discontinuities 
at e = —a, 0, and a in n(e) (see the insets in Fig. [2]) thus leading to A; = 2 solitons. A linear 
analysis around this state yields two unstable modes that exponentially grow with rates 271^2 in 
the first case and 27 it 2i// in the second. Which of the two cases is realized depends on the ratio 
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FIG. 2: (color online) fc = 2 normal soliton solutions to the Bogoliubov-de Gennes equations H]). Aq is 
the ground state BCS gap. At t ^ ±oo the 2-solitons tend to a normal state characterized by 2fc — 1 = 3 
jumps at e = and ±a in the average fermion occupation number n(e) (insets), a) is obtained from Eq. ([S]). 
where 71,2 are related to a and Ao by Eq. ([55]) . and 4.25a — Aq. b) corresponds to Eq. ([7]) with 4^ = Aq, 
47 = -^/iGa^ — Aq (see the text below Eq. ([5S])), and a — 3.75Ao. 01,2 = in both cases. For large 
separation, \a\ — 02] ^ 1, the 2-soliton solutions split into two single solitons (dashed lines), see Eqs. (|4]), 
(fT2| . (|90p. and Fig. [1] For |ai — q;2| -^ the two solitons merge into a single peak. In b) the amplitude of 
this peak is modulated with a frequency w ~ 4/^ = Aq as the two terms in Eq. ([7|) "rotate" with respect to 
each other. 

a/Ao- The values of 71,2 or 7ibz/i are also fixed by a and Aq, while 01,2 and ^\^2 are arbitrary real 
parameters. They control the separation between solitons and their relative phase in the 2-soliton 
solution, respectively. For sufficiently large \a\ — 02! Eq. ([5]) always breaks up into a sum of two 
single solitons of the form ([4]). The graph of | A(t)| in this case displays two well separated peaks, 
each representing a single soliton, see Fig. [2j 

The general fe-soliton solution of the normal type has the form (see Sec. IIV A|) 



|A(t)| 



-,2fc-l 



where D^ is the following determinant: 



Dr 
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Dk 



f(r~l) 



(8) 



(9) 



fir-l) ___ J2(r-1) 

j(m) jg ^YiQ rnth derivative of the function f{t) with respect to t, and 

^(cj)e-2»c,t 



2fc 



/(*) = E 



^n 



niy^j \'^j 



(10) 




FIG. 3: (color online) 3-normal-soliton solution to the Bogoliubov-de Gennes equations (H]) obtained from 
Eq. ([8]) with random Cj and (/>j. Aq is the ground state BCS gap. a) and b) differ only in the values of aj. 
In a) a2 — ai = 1.00 and as — a2 = 1.25; in b) a2 ~ ai — 10 and aa — 0:2 = 12. We see that for large 
differences between aj in b) the 3-soliton breaks up into a sum of three well separated individual solitons, 
see Eq. (|12p . while in a) the same solution but with small differences describes a complicated interference 
between the three solitons. 



The set of 2k complex parameters Cm ("frequencies" of the unstable modes) is complex conjugate 
to itself. Let us order this set so that c^+i = c*^ and Im(q) > for / = 1, . . . , A;. The constants 
A{ci) and A{ck-\.i) are related as follows 



A{ci) 



^ai+i<j}i 



A{ck+l) 



,-ai+i<f>i 



(11) 



where ai and (/)i are arbitrary real parameters. The single soliton ^ is obtained from Eq. ([5]) by 
setting k = 1 and ci = fj, + i^. The 2-soliton corresponds to A; = 2 and ci^2 = "^71,2 or ci^2 = i'J ± IJ^ 
to get Eq. ([6|) or ([7]), respectively. See also Fig. [3] for examples of 3-normal-solitons. 

At t ^ ±00 the /c-normal-soliton tends to a normal eigenstate that has at least 2A; — 1 discontinu- 
ities in the distribution function n(e). Linearizing the Bogoliubov-de Gennes equations around this 
state at large negative t, one obtains k unstable normal modes that grow as e"'^*'^'* for I = 1, . . . ,k. 
When the differences between parameters ai are large the /c-soliton solution ([5]) breaks up into a 
sum of k single solitons, 

k 

|A('=)(t,{Q,a„(/<J)| « J]|A«(t,Im(Q),a.)l, (12) 

where A^^'{t) denotes the A;-normal-soliton ([SD and |A''^)(t,Im(cj), aj)| stands for the single soliton 
^ with 7 = Im(ci) and a = ai. In this case, the plot of |A(t)| shows k well separated peaks 
(individual solitons) as illustrated in Figs. [2] and [3l Im(2Q) set the amplitudes and widths of 




FIG. 4: (color online) k = I anomalous soliton solution to the Bogoliubov-de Gennes equations ^ as given 
by Eq. ([13]). At t — > ±cx) the system tends to an eigenstate of the BCS Hamiltonian ([3]) with order parameter 
Aa- This eigenstate is characterized by 2/c = 2 jumps at e = ±a in the average fermion occupation number 
n{e) (inset). 7 and Aa in Eq. P^ are related to a and the ground state gap Aq via Eqs. (|15)) and (IST 
Here a = 1.47Aa and A^ = .09Ao. 



individual solitons, Re(2c/) are the frequencies with which they "rotate" with respect to one another 
as in Eq. ([7]), where Re(ci_2) = i/U, and ai and (j)i determine the separation between the solitons 
and their relative phases, respectively. 

B. Anomalous solitons 



Next, let us summarize the results of Sec. |V] for anomalous solitons. A single anomalous soliton 
has the following form (see also Fig. U]): 

2(f - Al) 



A(t) - A, 



Aa lb 7 cosh(At + a) 



(13) 



where A = 2y^7^ — A^ and a is an arbitrary real parameter as before. As t — > ±00 the state of 
the system tends to an anomalous eigenstate with the value of the BCS order parameter equal 
to Aa. In this state, all pairs in a certain energy interval around the Fermi level are excited, i.e. 



Em = y/f-m + Afl for \em\ < «• As a result, the distribution function n(e) has two jumps at e = ±a 
(inset in Fig. [3|). A linear analysis around this anomalous eigenstate shows a single unstable mode 
that grows as e . In general, 2k jumps in the distribution function of a stationary anomalous state 
give rise to up to k anomalous solitons. Note also that the anomalous soliton (J13p generalizes the 
normal one @ and turns into it when A^ = 0. 

A general /c-anomalous-soliton can also be constructed within our approach. However, here we 
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FIG. 5: (color online) k = 2 anomalous soliton solutions of the Bogoliubov-de Gcnnes equations ([T|), see 
Eqs. IM]) and (jlSp. a) corresponds to the choice of signs ++ and b) to — + in Eq. (|15p . At i ^ ±cx) 
both 2-anomalous-solitons tend to the same stationary state of the mean-field BCS Hamiltonian. In this 
state, n(e) has 2fc = 4 discontinuities at e = ±a and ±b (insets) and the anomalous average is equal to A^. 
Parameters 71^2 in Eq. (|15p and the value of A^ are determined by a, 6, and the ground state gap Aq, see 
Eqs. dSI]) and ^. For the above plots we used a = .87Aa, b = 4.27Aa, and IGA^ = Aq. 



present only an example of a 2-anonialous-soliton solution 

h{t) 



A(t) - A, 



h{t)h{t) - h'^it) 



(14) 



where 



hit) 



Aa , 71 cosh(Ait + ai) 72 cosh(A2t + 02) 



A]_A2 '^i('^i-'^2) ^2(^2 "-^i) 



(15) 



The lb signs can be chosen independently of each other. As in the case of the 1-anomalous-soliton, 
at large times the wave function asymptotes to an anomalous stationary state with order parameter 
Aa- This state has two unstable modes that grow exponentially with rates Ai^2 



The graph of the 2-anomalous-soliton solution ()14p consists of two peaks or dips depending on 
the choice of signs in Eq. (|15p . see Fig. [5j The parameters qi^2 take arbitrary real values. Their 
difference, \ai — 02!, determines the separation between the peaks in time. Similarly to Eq. ()12p . at 
large separations the fe-anomalous-soliton turns into a sum of individual solitons of the form (jlSh 

k 

(16) 



A^'\t,{K,ai}) -Aa^Yl [A(')(t, A,,a,) - A, 



i=l 



where A^^'[t, Xi,ai) is the single anomalous soliton ()13p with A = Aj and a = ai. 

The rest of the paper is organized as follows: in the next section, we review the basic setup of 
the problem and the tools (Lax vector and separation variables) necessary for deriving solitons. In 



10 



Sec. IIVI we perform linear analysis of equations of motion around normal and anomalous stationary 
states. This section also provides examples of normal and anomalous eigenstates that give rise to 
one and two normal and anomalous solitons. Sections IIVI and IVl are devoted to a detailed derivation 
of soliton solutions and a discussion of their main properties. 

II. REVIEW OF THE BASIC SETUP AND RELEVANT PREVIOUS RESULTS 

In this section we discuss the basic setup of the problem and introduce our notation (see Refs 



29 



30, l33|, M of the 



3d for more details). We also review the properties of the exact solution 

equations of motion needed for obtaining and analyzing the multi-soliton solutions summarized in 

the previous section. 

A. Notations and basic equations 

Here we review the model Hamiltonian and its mean- field equations of motion ([T|). The latter 
can be reformulated as equations of motion for classical spins (angular momenta) - this is the form 
we will be primarily using. We also describe normal and anomalous stationary states in terms of 
the spin variables. 



The dissipation 



Hamiltonian 2J, 



25 



ess dynamics of a fermionic superfluid can be modeled by the reduced BCS 



27 



38| 



^ = Y1 ^J'^i ~ ^ ^ cJtSV^J-^^T' (17) 

3 j,k 



st 



where hj = ^^.^i + Cj^Cjcr- This description is valid in the weak coupling regime at times shorter 
than the energy relaxation time r^ and for a system of size L smaller than the BCS coherence 
length ^. Under these conditions, the BCS order parameter is uniform in space, the interaction 
matrix elements can be evaluated at the Fermi energy ep yielding a single coupling constant g that 
is independent of j and k. The summations in Eq. (J17p over j and k are restricted to single particle 
energies \ej\ < D and \€k\ < D, where D is an ultraviolet cutoff for the pairing interaction. In 
metallic superconductors D ~ lod, where ujd is the Debye frequency. For atomic fermions D c^ ep. 
The offdiagonal interactions - terms of the form Cj^&j^^ c^iCr} with I ^ j or r ^ k can be neglected, 
since they are relevant only at times t > t^. 

The validity of the mean-field approach is rooted in the fact that each pair creation operator 
Cji^c-, in Eq. (fTTj) interacts with the collective pairing field g J2j CkiCki, which is expected to deviate 

11 



little from its quantum mechanical average A(t). For example, the mean- field is known to be exact 
'or the description of the low-energy properties of the Hamiltonian (|17p in the limit 5/Ao — > 



40, |45|, |47l], where 6 = (em+i — Cm) and Aq are the mean spacing between the single particle 
levels em and the ground state gap, respectively. Note that the conditions J <C Aq and L <^ S, are 
compatible in the weak coupling regime. 

We are interested in solving the Heisenberg equations of motion for Hamiltonian (J17p to de- 
termine the evolution of various correlators, e.g. {nm{t)), {cmi{t)cm'\{t)), and (c]„t(i)c|„i (0)- In 
mean-field approach, we replace the operator g'^jCki{t)ck^{t) in the Heisenberg equations with 
its quantum mechanical average 

A{t)=gY,{cki{t)ck^{t)). (18) 



Further, introducing 



we obtain 



Q 



2s^ = {flm) - 1, 



^in — ^m '^^rn — {CmlCmlfi 



(19) 



where A^, and —Ay are the real and imaginary parts of A = gYlim^m- In terms of Fm{t) = 
2isf^ and Grn(t) = 2is^, Green's functions at coinciding times, Eqs. (j20p are well-known Gorkov 



equations 



24 



4g]. The above procedure leading to Gorkov equations is essentially equivalent to 
taking the time-dependent wave function of the system to have the product form ([2]) at all times. 
Then, the Schrodinger equation takes the form of the Bogoliubov-de Gennes equations ([T|), which 



are in turn equivalent to Eq. (j20p with 

9„z _ IT/ |2 _ irr |2 - _ tt y* (9^^ 

Equations of motion (I20p are Hamilton's equations for the following classical spin (interacting 
angular momentum) model: 

n n 

H,i = Y,'^e^s]-9Y.'Pk^ (22) 

j=l j,k=l 



with the usual angular momentum Poisson brackets {s^, s^} = — s| etc. The summations in Eq. (j22p 
are restricted to the subspace of unblocked (with occupation numbers rij = 0, 2) single fermion 



12 



levels €j. The Hamiltonian (J17p does not have matrix elements connecting the unblocked levels 
to blocked ones {uj = 1). The latter are decoupled and their occupation numbers are conserved 
by the evolution. Note also that Eq. ([1]) conserves the norm |f/mP + iKnP = 1 and therefore the 
length of the spins is fixed \sm\ = 1/2. 

The normal and anomalous eigenstates of the BCS Hamiltonian discussed in the Introduction 
correspond to equilibrium spin configurations where each spin s^ is either parallel or antiparallel 



to the field br. 



40[. According to Eq. (j2ip . every such arrangement of spins uniquely determines 



an eigenstate of the form ([2]) and vice versa. The anomalous eigenstates yield 

where the x axis has been chosen so that the stationary value of the order parameter, A^, is real. 
The factor em = —1 if the spin is parallel to the field and Cm = 1 otherwise. The self-consistency 
condition Aa = g J^m ^m reads 

y /"^ =1. (24) 

This is the BCS gap equation, which determines the value of A^ in the anomalous state. The 
configuration of spins with all em = 1 is equivalent to the BCS ground state. In this case Aq = Aq 
- the ground state gap - and Eq. (j24p becomes in the continuum limit 



l"" ;/\, = | x = gi^FV^'-, (25) 

where vp is the density of states at the Fermi level and V is the volume of the system. In Eq. (j25p 
and throughout this paper we assume the weak coupling regime Aq <C D and a constant density 
of ej, z^(e) = I'F, in the continuum limit. A non-constant density of states modifies the value of Aq 
determined from Eq. (j25p but will not affect any other equations derived in the rest of the paper. 
As we will see, these equations are confined to energies of order Aq, while the density of states 
varies on an energy scale of order D or larger. Using Aq <^ D, we obtain from Eq. (I25|) 

Ao = 2De-^/^. (26) 

The configuration with only one flipped spin, e^ = — 1 and em^fc = l, corresponds to an excited 



state - it contains an excited pair and has energy 2w e| + Aq relative to the ground state. Similarly, 
having two spins parallel to the field is equivalent to an eigenstate with two excited pairs etc. 
Normal eigenstates are spin arrangements where each spin is along z axis, i.e. 

2s^ = ±1 = Im, s- = 0. (27) 

13 



They are also equilibria of the classical Hamiltonian ()22p according to Eq. (j20p . The Fermi ground 
state corresponds to Im = — sgn em, while other normal eigenstates can be obtained from this state 
by moving pairs of fermions from levels below the Fermi energy to levels above it and flipping the 
corresponding spins. 

B. General properties of the dynamics 



In this subsection we introduce the Lax vector construction 30|, |33|] , which plays a central role 
in analyzing the dynamics of the BCS Hamiltonian. We also define the separation variables and 
describe the general features of the dynamics. 

The dynamics of the classical Hamiltonian (j22p or, equivalently, Eqs. ()20p and ([T|) turn out to 
be integrable. A convenient tool for their analysis is the Lax vector defined as 

n 
m=l 

where u is a complex parameter, z is a unit vector along z axis, and n is the total number of spins. 
The length of this vector is conserved by Eqs. (j20p for any u, i.e. 



dt 



0. (29) 



M] for Eqs. ([20]). For 



For this reason L^(ii) can be viewed as the generator of the integrals of motion 

example, its zeroes are conserved and constitute a set of independent integrals. Another possible 

choice for the integrals is e.g. the set of the residues of L^(u) at the poles at u = em- Note that 

where Q2n{u) is a (spectral) polynomial of order 2n. We also have 

L2(n) = Ll{u) + L^{u)L+{u), (31) 

where L-{u) = Lx{u) — iLy[u) and L^^y^z are the components of the Lax vector L(ti). 



48|, 



49(] in which Eqs. ([20 



To obtain solitons, we need to introduce new dynamical variables u„ 
separate and can be integrated. The separation variables are defined in terms of the "old" dynam- 
ical variables Sj as solutions of the following equation: 



L.{Um)=Y^ ^ = 0- (32) 
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This equation has n — 1 solutions since, when L-{u) is brought to a common denominator, its 
numerator is a polynomial of order n — 1. Consequently, there are n — 1 separation variables Um- 
Eq. ()32p can be inverted to obtain the spins in terms of the separation variables as follows 

where J^ = Jx — iJy as usual and J = J2j ^j ^^ ^^^ total classical spin. 
In terms of the separation variables Eqs. ([20}) read 



«7 = ^-tt^^t^^' (^^) 





'^^^/Q2n{uj) 


1 n 




rim^iK-^m)' 




,L 


/ n 





1, (34) 

(35) 

. j=l m=l 

An important observation [34,] is that main properties of the dynamics can be effectively dis- 
cerned by analyzing the zeros of L^ {u) . According to Eq. (|30p , these are the roots of the spectral 
polynomial Q2niu) and we will often refer to their configuration in the complex plane as to the root 
diagram of L^(ti). Since Q2n{u) is positively defined, it has n pairs of complex conjugate roots. 
For generic initial conditions all 2n roots are distinct. In this case the dynamics of the system 
is quasiperiodic with n incommensurate frequencies and any dynamical quantity, e.g. the order 
parameter A(t) = gj^{t), typically contains all n frequencies. 



Significant simplifications occur when some roots are degenerate [29l. l33l| . It is important to dis- 
tinguish between real and complex double roots. Note that any real root of Q2n{u) is automatically 

a double root (zero) because Q2n{u) is positively defined. A real zero c of L^(ti) must also be a 

i I 

zero of all three components of L(ti)[50]. Further, note from Eq. ()32p that one of the separation 
variables must coincide with c. In other words, it must be time-independent as it is "frozen" into 
the real root c. Eq. (j34p shows that this is an allowed solution of the equations of motion for 
the separation variables. This freezing of a separation variable can be translated into a genuine 
reduction of the number of degrees of freedom by one so that the dynamics of the Hamiltonian ()22p 
with n spins reduces to that of the same Hamiltonian but with n — 1 spins. In general, n — m real 
zeros (or equivalently 2m complex zeros) mean a reduction of the dynamics to that of 2m effective 



spins, see Ref. 



29|and 



331 for details. Below we will often encounter a situation when L^(u) has a 
number of real zeros and consequently a number of separation variables are frozen. The remaining 
variables we call unfrozen. 

Let Un-i = c be the separation variable frozen into the double zero of Q2n{u). Consider Eq. (f34l) 
for j 7^ n— 1. Both the numerator and the denominator of the right hand side contain a factor u — c, 
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which cancels lowermg the order of the polynomial under the square root by two. Suppose Q2n{u) 
has n — 2k double real zeros. Then, there are 2k — 1 unfrozen separation variables ui, . . . , U2k-i- 
For these variables Eq. ([M|) can be brought to the following form^] with the help of Eq. ([30]) : 

2'="^ u^du- 

Y, ' ' p= = ^i9dt6i,2k-2, 1 = 0,..., 2k -2, (36) 

where L^(ti) is obtained from L'^(u) by removing all real zeros Cm, i.e. 

L [u) 



2^ 

va) 



nm(^-c« 

III. LINEAR ANALYSIS AROUND STATIONARY STATES 

In this section, we analyze equations of motion linearized in the vicinity of normal and anomalous 
stationary states. We show that the separation variables Uj are the normal modes of the linearized 
problem. Some of the stationary states are unstable. As we will see in the next section, the 
corresponding normal modes become solitons in the nonlinear regime. 

A. Frequencies of oscillations around stationary states 

Here we show that the frequencies of small oscillations around normal and anomalous states are 



determined by the zeros of L (u), see also Ref. 



34l . When one of the frequencies is complex, the 



state is unstable and the corresponding mode grows exponentially. 

The linear analysis of equations of motion (j20p around stationary states greatly simplifies in 
terms of separation variables. According to Eq. (j34p . stationary positions of Uj are the roots of the 
polynomial Q2n{u) (or equivalently the zeros of L^(u), see Eq. (|30|) ). Let us determine the form 
of L^(u) in the stationary states. Consider first the anomalous states (j23|) . Using Eqs. ([23]) . (p8]) . 
and ()24p . we obtain 

L(n) = (A,x - uz)Ls{u), L\u) = {u^ + l^l)Ll{u), (37) 

where x is a unit vector along x axis and 

n 

Ls{u) = y """^^ e^ = ±1. (38) 



Note that when the right hand side of Eq. (j38p is brought to a common denominator, the numerator 
is a polynomial of order n — 1. Therefore, L^(u) and consequently Q2n{'u) have n — 1 double zeros 

16 



Cr — the solutions of the equation Ls{cr) = 0. In addition, we see from Eq. (j37p that there are two 
roots u = ±iAa, i.e. 

n-l 
Q2n{u) = {u" + A2) J{{u- Crf. (39) 

r=l 

When the spins Sj deviate from their equihbrium positions (j23p . the roots c^ of polynomial 



Q2n{u) shift to Cr + 6cr 5l|] . Linearizing Eq. (j34p in deviations (Jc^ = a^ + i&r and (^Ur = n,- — c,. — a^ 



around the stationary positions Ur = Cr and using Eq. (j39p . we obtain 

<5u, = 2iVc2 + A2^(<5u,)2 + 62 (40) 



with a solution (Ju^ = ^r sin[(j-'r'(i — ^o)]) where ujr = '^\/c'^ + A^. Linearizing Eq. (j33p . one derives 
the spin variables in terms of 5ur. At this point we are interested only in the frequencies uJr- 

We conclude that the separation variables are indeed the normal modes of the linearized problem 
(since they contain a single frequency). The frequencies of small oscillations around anomalous 
stationary states are related to the double zeros c^ of ^"^{u) as uJr = 'i'^J c^ + A^. If any of Wr has 
an imaginary part, the stationary state is unstable. 

Next, consider linear analysis around normal eigenstates (j27p . In this case all spins are along z 
axis. It follows from Eqs. (I27p and ()28p that L(n) = Ln{u)z, where 

1 " /• 

We see that all zeros Cr of L^(u) = L'^{u) are double zeros. There are n of them as the numerator of 
Ln{u) is a polynomial of order n, i.e. Q2n{u) = Y[r=ii'^ ~ '^r)'^- ^^ before, the stationary positions 
of separation variables are Ur = Cr- Note however that there are only n — 1 separation variables, 
so one of the n zeros Cr must remain vacant. 

Next, we show that the frequencies of small oscillations around a normal stationary state are 
iOr = 2cr- When one of the zeros Cj- is complex, the oscillatory behavior is replaced with an 
exponential growth, i.e. the stationary state is unstable. Note that for small deviations from a 
normal state the xy components of the total spin J are small. Therefore, in linear approximation 
we can set the separation variables Uj to their equilibrium values in Eqs. (f33]l and (f35]l . Uj = Cj, 
i.e. only J_ is time-dependent. As mentioned above, Ln{u) has a vacant zero (say Cr) which does 
not correspond to any separation variable. Eq. (|35p yields 



d(ln J- 
2idt 



n J n 

j=l m=l 

17 



+ Cr = Cr, (42) 



where we used the fact that the contribution in square brackets vanishes. This can be seen by 
observing that, since Cm are the zeros of Ln{u), Eq. (|4ip can be written as 

Ln{u) = ~-^f^^^^. (43) 

Expanding the right hand sides of Eqs. (j43p and (j4ip in 1/n, matching the coefficients at 1/u, and 
using 2Jz = Ylj h (this follows from Eq. (I27p ). we see that the sum of terms in square brackets in 
Eq. (j42p is indeed zero. It follows that J_ oc e^^*^''* and from Eq. (ISSp we also derive s^ oc e"^*'^''*. 
Thus, the frequencies of oscillations around normal stationary states are Ur = 2cr. 

B. Examples of root diagrams of L^(u) 

Here we provide examples of root diagrams - configurations of solutions of the equation L-^ (u) = 



in the plane of complex u, see Ref. |34 for more examples. We saw that the zeros of L'^(u) 
evaluated in stationary states determine the frequencies of oscillations around them. Moreover, 
the most important features of the dynamics, e.g. the behavior of A(t) at large times, can be 
predicted by inspecting the root diagramfsj], see also the discussion below Eq. (f35l) . Similarly, we 
will see that the root diagram determines the number and properties of solitons corresponding to 
a given stationary state. 

For simplicity, we assume particle-hole symmetry, i.e. the single fermion energies {em} are 
symmetric with respect to zero (Fermi level). According to Eq. (I23p . this means 

(44) 

where s^ = s(em)- These relations can also be derived from Eq. (J19p using particle-hole trans- 
formation for fermion creation and annihilation operators Ca-{—em) '^ Ca{em)- Note that relations 
([44ll are preserved by equations of motion (f20j) and also imply Ay{t) = 0. 

1. BCS ground state 

As discussed below Eq. (pi]) , the BCS ground state corresponds to e^ = 1- We note from 
Eq. ([23]) that spin components s^{em) and s^(em) in this state are continuous functions of single 
particle energy em- It follows from Eqs. (f37|) and ([38]) that L^(u) = has two single roots at 
u = ±iAo, see Fig. [H and n — 1 double roots Ck that are the solutions of the following equation: 

n 

LJu) = y , = 0. (45) 

^,2(n-6m)V4:TA^ 
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FIG. 6: (color online) The roots of L^(u) = (the root diagram) for the BCS ground state in the complex 
u plane. There is a line of double real roots (zeros of L^(u)) from —D to D, where D is the high-energy 
cutoff on the single- fermion states participating in the BCS Haniiltonian ()17|1 . In addition, there are two 
imaginary single zeros ±iAo, where Aq is the ground state gap. Frequencies lu of small oscillations around 



the ground state are related to the zeros c a,s u; — ^ (? + Aq, see Sec. IIIIB 11 We have a;(e) = i/e^ + Ag, 
where — Z? < e < D, and uj = 0. The inset shows the spin component s^(e) in the ground state. Since it has 
no discontinuities, there are no complex double zeros. 



All n—1 solutions are real. This can be seen by noting that Ls{u) changes sign between consecutive 
em- Indeed, let e^ be ordered so that ei < 62 < • • • < e„. Since Ls{u) — > -|-oo as u — > e+ and 
Ls{u) — > — cxD as u ^ ^m+i' there is a point u = Cm in the interval (em,em+i) where Ls{cm) = 0. 



According to the previous subsection, this yields a frequency ujm = 2y^c^ -|- Ag of small oscillations 
around the BCS ground state. In the continuum limit, when level spacings em+i — Cm tend to zero, 
we have Cm ~ Cm and u>m ~ '^\l ^ + ^o- ^"^ this limit, Cm densely fill the interval from — D to 
D = max I Em I, i-e. L^('u) has a line of double roots as shown in Fig. [6l Frequencies Wm are also 
the energies of excited pairs - excitations obtained by flipping the spin s^ from its ground state 
position antiparallel to the held b^ = (— 2Ao, 0, em) to an equilibrium position parallel to the held 



bm, see Ref. 
spectrum. 



4d for a discussion of this relationship between the frequencies and the excitation 



2. Excited anomalous states 

Next, consider two examples of anomalous stationary states obtained from the BCS ground 
state by flipping spins in certain energy intervals. 
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FIG. 7: (color online) Zeros of L^(?i) for an excited anomalous state in the complex u plane. This anomalous 
state has 2fc — 2 jumps at e = ±a in the spin component s^(e) (inset) leading to 2k double imaginary zeros 
±17. In addition, L^(ii) has a line of real double zeros from —D to D and single zeros ±iAa, where Aq is 
the value of the order parameter in this state. In the present case, there is a single (k — 1) unstable mode 



that grows with the rate a/7^ — A^ giving rise to a single anomalous soliton shown in Fig.lH see Sec. lIIIB^ 



Example 1 

First, let spins in the interval {—a, a) be flipped, i.e. e^ = sgn(|em,| — cl). This means that 
the Cooper pairs in this energy interval are excited [38l.l4Cll|. Eq. (j23p implies that spin components 
s^(e) and s^{e) are discontinuous at e = ±a (see the inset in FiglT]). As before L^(n) has two 
single zeros at ti = iiA^ and n — 1 double zeros Ck that are the solutions of Ls{ck) = 0, where 
Ls{u) is given by Eq. ([38|) . The difference is that in this case two of Ck can be imaginary. Suppose 
Cmi < —a < emi+i and 6^2 < o < em2+i- Then, e^i = 1 while Cmi+i = —1 and similarly for 1712 
and we are no longer guaranteed real zeros of Ls{u) in intervals {emi, e-mi+i) and (6^2; ^m2+i) as in 
the BCS ground state. Instead, Ls{u) can acquire two complex conjugate zeros. In the particle-hole 
symmetric case Ls{—u) = —Ls{u), which implies that these zeros must be purely imaginary as in 
Fig. [71 The remaining n — 3 double zeros of L-^ (n) are real and lie between consecutive ej . In the 
continuum limit, they merge into a continuous line of double zeros between —D and D as in the 
ground state. 

To determine the two imaginary zeros c = ±27 in the continuum limit, we rewrite the equation 
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Ls{u) = in the form 

sgn(e — a) de 



0, (46) 



/o (e2+72)y^2^A2 
where we used Ls{—u) = —Ls[u) and took the ultraviolet cutoff D to infinity. In terms of 



^, . I de 1 , 

F{e) = / ^ = ^ In 



7y^^TA2 + ev/^23A2 



Eq. (|16]) reads F(+oo) = 2F{a). This equation has a unique positive solution 



(47) 



7=la+\/^^TAl)— . (48) 



Note however that the gap equation (j24p has solutions only for sufficiently small a. To see this, 
we write down Eq. (I24|) for the order parameter A^ in the anomalous state where Cm = sgn(|em| —o) 
and for the gap Aq in the BCS ground state where e^ = 1- Equating the left hand sides of the 
two equations, we obtain 

"^ de _ f^sgn{e-a)de 



/o y?TA^ Jo \/?ta2 ■ 

In the D ^ oo limit this equation yields 

A^ - 2AoA2 + A^A, -4a^ = (50) 

together with the condition A„ < Aq. 

The analysis of Eq. ([50]) shows that there are two solutions Aq < Aq provided 3\/3o < Aq and 
no solutions otherwise. Interestingly, for one of the solutions 7 < A^, while for the other 7 > A^. 
Which solution do we choose? Note that the quantum Hamiltonian (J17p has 2" unblocked states. 
Correspondingly, there are 2" choices of e^ = ±1. More than one solution for a given selection 
of Cm means that we have more states in the mean-field than there are eigenstates of the original 
quantum Hamiltonian. It is natural to expect that among the two solutions for A^ the one that 
yields a stable anomalous state corresponds to the quantum eigenstate. We have shown above 
that frequencies of small oscillations around anomalous stationary states are related to the zeros 



Ck as u)k = \J c\ + A^. For the zeros ±Z7 we have oj^y = iy^7^ — A^. We see that for 7 > Aq the 
frequency is imaginary and the corresponding normal mode grows exponentially. Therefore, the 
solution Afl < 7 yields an unstable anomalous state, while for A^ > 7 we get a stable state. Both 
states however can play an important role in the description of the dynamical problem ([1]) . 
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FIG. 8: (color online) Zeros of L^(?i) for an excited anomalous state in the complex u plane. This anomalous 
state has 2k — 4 discontinuities at e = ±a and ±6 in the spin component s^(e) (inset) leading to 2fc double 
imaginary zeros zLiji,2- In addition, L^(u) has a line of real double zeros from —D to D and single zeros 
±iAa, where Aa is the value of the order parameter in this state. In the present case, there are k = 2 



unstable modes that grow with rates \/jf 2 ~ A^ giving rise to 2-anomalous-solitons shown in Fig. [5l see 
Sec. lIIIB2l 

Example 2 



A more involved example of an anomalous state is obtained by flipping spins in two energy 
intervals, e.g. in intervals (—6, —a) and (a, b) symmetric with respect to the Fermi level. This 
implies em = sgn(|em| — a)(|em| — b). Now the spin components have fom' discontinuities at e = ±0 
and e = ±6 (inset in Fig. [8l). Correspondingly, L^(u) can have four complex double zeros ±171^2 
as in Fig. [8] in addition to two single zeros at u = iiA^ and n — 5 real zeros on the line from —D 
to D. This follows in a manner similar to the above analysis of the state with Cm = sgn(|em| — «)• 
In general, 2k discontinuities in spin components in an anomalous stationary state can lead to 2k 
complex double roots. 

To determine the complex zeros in the continuum limit, we repeat the procedure that lead to 
Eq. (j48p . Now we derive 2F{b) — 2F{a) = F{+oo), where F{e) is given by expression (|47p . This 
equation has solutions ±^71 ^2) where 



Aa 



xy 



l±.J{xy + lf-A{x + y-l), 



-A2 



+ A2 



and yA^ 



6^ + A^ + b\ . The gap equation 



(51) 

in terms of x and 
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y takes the form 



xy 



Ao 

A-- 



(52) 



where Aq is the ground state gap. Using these equations, it is not difficult to select a and b so that 
7i 2 are real and 72 > 7i > A^ as shown in Fig. [8l This is the choice we will need in Sec. IV B[ 
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FIG. 9: (color online) Zeros of L^(w) for the Fermi ground state in the complex u plane. The single 
discontinuity (fc = 1) in the spin component s^(e) at e = 0, see the inset, leads to a single (2fc — 1 = 1) pair 
of complex double zeros at ±jAo/2, where Aq is the order parameter in the BCS ground state. There is also 
a line of double real zeros stretching from —DtoD, where D is the high-energy cutoff on the single- fermion 
states participating in the BCS Hamiltonian (|17p . Frequencies lu of small oscillations around this state are 
related to the zeros c as u; — 2c, see Sec. lIII A1 In the present case, there is a single {k = 1) unstable mode 
that grows with the rate 7 = Aq giving rise to fc = 1 normal-soliton shown in Fig. [T] 



3. Fermi ground state 

We saw that in normal stationary states all zeros of L^(u) are double degenerate and are 
solutions of the equation Ln{u) = 0, see Eq. (j41|) . The Fermi ground state has all states below 
the Fermi energy occupied and states above it empty. This corresponds to 2s^ = Ij = — sgne^. 
Therefore, the zeros are determined by the following equation: 






sgn €j 

U — Ej 



2 
9' 



(53) 



There are n solutions each one being a double zero of L^ {u) . The analysis of Eq. (|53p is similar to 
that of Eq. (|15|) . Eq. ([53l) has real roots between consecutive ej except when sgn e^ changes from 1 
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to —1. Therefore, there is a real root Cj in each interval (e^, Cj+i) except for the interval containing 
the Fermi level. Since there are n — 2 such intervals, n — 2 roots are real while the remaining 
two can be complex. Due to the particle-hole symmetry (j44p the complex roots must be purely 
imaginary, see Fig. [9l They also must be complex conjugate to each other as Eq. (j53p is invariant 
under complex conjugation. 

In the continuum limit the spacing between ej vanishes and for the real roots we have Cj ~ ej, 
i.e. L^(u) has a line of double real zeros stretching from —D to D, Fig. [9j To determine the two 
imaginary roots ±17, we rewrite Eq. (j53p in the integral form 

Ve 



de 



+ 



de 



2 
A' 



Z7 € + IJ 

Using Eq. (126]) . we obtain in the weak coupling regime Aq <^ D 



7 



Ao 
2 ' 



(54) 



Thus, according to the discussion in the previous subsection, equations of motion (j20p linearized 
around the Fermi ground state show n — 1 stable modes with oscillation frequencies ujj ~ 2ej and 



pAot 



u 



4l[. Note also that the z component of spins 



one unstable mode that grows as e^"^* 

2s^{ej) = —sgncj in the Fermi ground state experiences a single jump at the Fermi level. 
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FIG. 10: (color online) Zeros of L^(u) for an excited normal state characterized by 2fc — 1 = 3 jumps in the 
spin component s^(e) aX e — ±a (insets). There are k = 2 pairs of complex conjugate double zeros (identified 
with crosses) leading to fc = 2 normal-soliton solutions shown in Fig. [2l The complex zeros are determined 
by a and the ground state gap Aq, see Eq. (|55|). a) corresponds to the case 4a < Aq and Fig. [2^) and b) 
corresponds to 4a > Aq and Fig.I^b). 
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4- Excited normal state 

Now consider a normal stationary state where s^(ej) has three discontinuities. We require 
2s^{ej) = —1 (1) for large positive (negative) €j. Otherwise, the first term in Eq. (j22p is not 
minimized at large Cj and single particle states far from the Fermi level are affected by the pairing 
interaction, which is unphysical. Under these conditions the total number of discontinuities in 
s^{ej) must be odd. Therefore, the next option after the Fermi ground state that has one jump is 
a state with three jumps in s^(ej). 

Let 2sj = Ij = — sgnej(e^ — a^), i.e. spins in the interval \ej\ < a point in directions opposite 
to those in the Fermi ground state, see the insets in Fig. [TOj The solutions of L^(m) = are 
determined in the same way as for the Fermi ground state. In the present case, we find that there 
are n — 4 double roots located between €j and e^+i except when Ij and /j+i have different signs. 
The remaining four double roots can take complex values. In the continuum limit, L^(u) has a line 
of double zeros from —D to D and four isolated complex zeros c = ^71, 2 and c = —^71, 2 shown in 
Fig. [TOl where 



Correspondingly, there are two unstable modes (one for each pair of complex conjugate zeros). If 
a < Ao/4, Eq. (j55p yields real 71^2 (Fig. llOb) and the unstable modes grow as e^"^^* and e^^^*. For 



a > Ao/4 we have c = ±/x it ij (Fig. [TOb). where /x = Ao/4 and 7 = y^a^ — Aq/16. In this case 
unstable modes diverge in an oscillatory manner as e^^*'^*e^'^*. In general, a normal stationary 
state with 2k — 1 discontinuities in s^{ej) is characterized by up to 2k complex double zeros of 
L^(n) and k unstable modes. 

IV. NORMAL SOLITONS 

In this section, we determine solutions of equations of motion ([T]) that asymptote to normal 
stationary states at t — > ±00. In particular, we derive equations ^ - (J12p for normal solitons, see 
also Figs. [H O and [3l That these solutions are solitons can be seen in a number of ways. First, 
these are trajectories that connect an unstable equilibrium to itself, i.e. they start in an unstable 
stationary state at t — > — 00 and return to it at t ^ cxd. This is typical of solitons [43], see the 
paragraph preceding Eq. dJ]). Second, as we will show, in a certain regime the solution splits into 



a sum of single solitons as it should 
are in terms of elementary functions 



■3\. Finally, in contrast to the general solution these solutions 

44|. 
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A. General fc-normal-soliton solution 

Consider a general normal stationary state with 2k — 1 discontinuities in s^(ej). Suppose 
L^(n) = L'^{u) has 2k complex double zeros ci,C2, ■ ■ ■ ,C2k, i-e. there are k unstable modes in 
the linear analysis. Let us solve equations of motion for separation variables (j36p for this state. 

First, we derive a useful equation for the time-dependent gap function A(t). Ln{u) has 2k 
complex conjugate zeros ci,C2, ■ ■ ■ ,C2k and n — 2k real zeros. Bringing Eq. (j4ip to a common 
denominator, we obtain 

r ( N 1 P2k{u)Rn~2k{u) , . 

Lniu) = .=— — , (56) 

where P2k{u) = Ylr=ii^ ~ ^r) and Rn-2k{u) represents the contribution of the real zeros. Both 
these polynomials have real coefficients. Eq. ([3TI1 yields 

[L,{u) - Lniu)] [L,{u) + Ln{u)] = -L^{u)L+{u) (57) 

This equation implies 

r ( -.^r t . 2 Sk{u)Sl{u)Rn-2k{u) 

L^(u)+Ln{u) 



9 Umi^ - ^m) 



. . ^ r ^ X 9 J J Tk-i{u)T*^^{u)Rn-2k{u) 

L:,[U) - Ln{u) = --J^J^ 



(58) 



2 " Umi^-e. 



r ^ ^ T Skiu)Tk-liu)Rn^2kiu) /^-ox 

^-"'' = -'- n,„("-^.,.) — ■ <'"' 

where Sk{u) and Tk-i{u) are polynomials in u of orders k and k — 1, respectively. The coefficient 
at highest power of u is equal to unity in all polynomials. The coefficients of 5'^('u) are complex 
conjugate to those of Sk{u) and similarly for T^_^{u). The prefactors in Eqs. (fSHIl and ([59ll are 
obtained from large u behavior. For example, Eqs. (j28p and (14ip imply Lz{u) + Ln{u) ~ —2/g at 
u — > oo. The right hand side of the first equation in (j58p has the same large u asymptote. The 
polynomial Rn-2k{u) is common to all components of L(ti) since any real zero of L^(n) is also a 
zero of Lx^y^z LSfl]. Subtracting the second equation in ([58]) from the first one and using Eq. ([56]) . 
we derive 

P2k{u) = Sk{u)Sl{u) + \^-Tk-x{u)n_^{u), (60) 

where we used A(t) = gj^{t). We will need this equation below to determine |A(t)|. 
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Now consider equations of motion (j36p . It follows from Eq. (j3ip and the definition (j32p of 
separation variables Uj that L^(uj) = L1{uj). Comparing the large u behavior of syT^iu) — > 1/g 



and Lz{u) -^ — 1/5, we conclude that y^L^ {uj ) = —Lz{uj). Eq. (J36ll takes the following form: 



2^-1 u'dn 



V .3 ^ — = 0, / = 0, ...,2A:-3, 

,=1 -^^(%)nmK-e'n) 



(61) 



2'=-! nf-^du 



V ^ ^ = -2igdt, 



where Lz{u) = Lz{u)/Rn~2k- Eqs. (f32]l and ([59]) imply that unfrozen separation variables 
ui, . . . ,U2k-i are the roots of Sk{u)Tk-^i{u). Let ui, . . . ,Uk-i be the roots of Tfc_i(n) and 
Uk, ■ ■ ■ ,U2k-i the roots of Sk{u). Eq. ([^Sj) reads Lz{uj) = Ln{uj) for j = 1,...,A; — 1 and 
Lz{uj) = —Ln{uj) ioT j = k, . . . ,2k — 1, where L„(u) = Ln{u)/Rn-2k- Further, using Eq. ([56]) . we 
obtain from Eq. ()6ip 



2A_^i u'dti- ^ -u'dn,- 

E7jVt-E^Vt = -2^'^*'^''2'^'2, / = 0,...,2fe-2. (62) 

f^l P2k{uj) j-lP2k{Uj) 

Eq. (j62p does not contain square roots in contrast to Eq. (j36p and can be integrated in elementary 
functions. To do so, we expand the ratios v!" / P2k{u) in elementary fractions 

irVr = y 1 TtT^^^i V' ^ < 2A:. (63) 

^2fc («) ^^ (n - c„) n^ym (cm - 9 ) 

This identity can be verified by comparing residues at poles n = c^ on both sides. Using expansion 

(j63l) in Eq. (I62|), we obtain 

2fc i , 

^^llljj^m\Cm Cj) 

where / = 0, . . . , 2/;; — 2 and 

. , Uj Cm . -, Uj Cffi -'fc— l(,Cmj -'-k—ly^mj 

j=k ■' 3=1 ■' 



Integration of Eq. (I64p results in 

2k 



En 7 ""_ . =-2i^^^,2;c-2 + i^(Q), Z = 0,...,2fe-2, (66) 

where E{ci) are the integration constants. These equations are linear in Xm with the general 
solution 

Xm = -2icmt + E{cm) + G{t). (67) 
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E{cm) are new time-independent constants and G{t) is an arbitrary function of t. 
Using the definition of Xm in Eq. (i65]l , we find 

JkyCm) 



-A{cm)F{t)e 



^ICfnt 



m = 1, . . . , 2k, 



(68) 



Tk-l{Cm) 

where A{cm) are complex constants and F{t) is a function of time to be determined below. Eqs. (j68p 
are 2k linear equations for 2k — 1 coefficients of polynomials Sk{u) and Tfc_i(n). The compatibility 
condition for this linear system yields a linear equation for the function F{t). We derive 



F{t) = (-1)'^2 



fcr,2fc-2 



Dk-1 
Dk 



(69) 



where the determinant D^ is given by 



Dr 



f 



... /(-I) 



(70) 



J(r-l) ^^^ J2(r-1) 

f^^' is the jth derivative of the function f[t) with respect to t, and 



2fc 



m = E 



A{cm)e 



■^^Cttt, t 



(71) 



fl n«^m(Cm - Q) ' 

To relate |-F(t)| to |A(t)|, we use Eq. (jGOp . This equation also imposes certain restrictions on 
complex constants A(cm)- Setting n = Cm in Eq. ()60p and using the fact that Cm are the roots of 
P2k{u), -P2fc(cm) = 0, we obtain 

T{Cm)T*{Cm) 4 ■ ^' ^ 

Note that while the coefficients of the polynomial S^{u) are complex conjugate to those of Smiu), 
S*{cm) is not complex conjugate to S{cm) since Cm is complex. Instead, we have S*{cm) = [Sic*^)]* , 
i.e. S*{cm) is conjugate to S{c'^). Using this and Eq. ([68]) . we obtain from Eq. ([72]) 

|A(t)|2 



A{c^)A*{0\F{t)\' 



(73) 



where j4*(c^) is the complex conjugate of A{c'^) - the constant corresponding to the zero cj^ 
complex conjugate to Cm (recall that the zeros Cj of L^(u) come in complex conjugate pairs). 
Eq. (f73|) implies that the product A{cm)A*{c'^) is independent of m. With no loss of generality we 
set 



A{cm)A*{c: 



-1. 



(74) 
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Any other real value will rescale |-F(i)| without affecting |A(t)|. Therefore, we have |A(t)| = 21^(^)1 
and 



|A(t)| = 2 



2k-l 



Dk- 



Dk 



It follows from Eq. ()74p that the constants A{cm) can be parameterized as follows 



A{ci)=e^'+'^\ A{ck+i) = -e 



-ai+i<f>i 



(75) 



(76) 



where ai and (j)i are arbitrary real parameters and we ordered the 2k zeros Cm so that Ck+i = c* 
and Im(q) > for / = 1, . . . ,k. 

Eqs. ([75]) . ([7n|) . ([7T]) . ([55]) and ([75]) fully describe the general /c-normal-soliton solution (examples 
for k = 1, 2, and 3 are shown in Figs. [H [21 andO respectively). They contain 2k zeros Cm fixed by 
the normal stationary state corresponding to this solution. This state has 2k— 1 discontinuities in 
s^(em)- The zeros Cm are the roots of the equation Ln{u) = 0, where Ln{u) is given by Eq. (J4T]) . 
The 2k real parameters ai and (pi in Eq. ()76p are arbitrary. That the general /c-normal-soliton 
should be indeed characterized by 2k arbitrary real parameters is seen from the discussion in the 
paragraph following Eq. 



). As mentioned there (see Refs. |29| and l33l for details), a real double 
zero of L'^(n) effectively reduces the number of degrees of freedom (spins) by one. Since in the 
present case we have n — 2k such roots, it can be described by 2k effective spins. Then, there are 
4A: initial conditions (two angles per each spin). 2k of these are determined by the 2k integrals of 
motion Cm, while the other 2k correspond to ai and (pi. 



B. Matching soliton constants to spin configuration at large negative time 



Here we show that the /c-soliton (j70p tends to a normal stationary state in f — > ±00 limits and 
relate the constants ai and (pi to the deviations of spins from this state at large negative times. 

First, let us evaluate expression (fTOJ) for large negative t. To this end, we keep in Eq. (fTTI) only 
the exponents that diverge in the t -^ —00 limit, i.e. the k terms that have Im(cm) < 0. After 
some manipulations with the rows of determinants Dk and Dk-i, we derive 



|A(t)| 



E 



„ — 2iCmt pam+i4>m ("^ 



^m) iliV^m Cj j 



(77) 



„=1 lli^mVCm Cj) 

We see that |A(f)| oc e"^'*'* at large negative i, where 7 is the minimum of |Im(cm)|. Quantities 
J±{t) and F{t) behave in the same way as they are proportional to |A(t)|. According to Eq. ()68p 
as t — > —00 either Sk{cm) ^ or Tk-i{cm) -^ except for the zero Cm with Im(cm) = ^7- Since the 
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unfrozen separation variables are the roots of either Sk{u) or Tk~i{u) (see Eqs. ([32]) and ([59]) 1. they 
must tend to their stationary state positions Cm- Observe also that since J±{t) — > 0, the second 
equation in (j57p and Eq. (|59p mean Lz{u) — > Ln{u) and L^{u) -^ 0. It follows from the definition 
(j28|) of L(n) and Eq. (j4ip that s^'^ -^ and Sy — > ^j/2. Thus, the fc-normal-soliton tends to the 
normal stationary state that has the same values of zeros q. The analysis of the t ^ oo limit is 
completely analogous. 

Next, consider the limiting stationary state. There are 2k zeros Cj and only 2k — 1 unfrozen 
separation variables, i.e. one of the zeros Cj (say Cr) remains vacant. Suppose spins deviate from 
this stationary state keeping the values of integrals of motion a the same. Since J_ = in normal 
states, Eq. ([33|) yields to the linear order in the deviation 

where ul = Ci are the stationary positions of the unfrozen separation variables and Rn~2k{^j) is 
the contribution of the frozen ones. The frozen variables are located in real zeros of L„(u), which 
are also the zeros of Rn-2k{u), see Eq. ([^U|) . Further, Eqs. (jlT|) and ([^H]) imply 

1 Y^ Ij _ I Rn-2k{u){u - Cr)Ui{u - uf^) ,„„. 

g^2^^2{u-e,) 9 UJu-e^) ' ^'^ 

Equating the residues at poles u = ej on both sides, we obtain 

Rn-2k (gj ) Ui (gj - uf^ ) _ lj9 

Substituting this into Eq. ([78|) . we find 



ij 



(80) 



Finally, J-(i) is determined from Eq. (j42p . which was also derived in a linear analysis around 
normal stationary states. The difference is that there we considered generic deviations when the 
integrals of motion Cj also deviate from their stationary state values. Nevertheless, Eq. ()42p is the 
same in both cases and integrating it, we obtain 

A{t) = gJ^it) = Pre-'^'^'-', (81) 

/ -2iCrt 

These are particular solutions of the linearized equations of motion. They describe an unstable 
mode with complex "frequency" 2cr. 
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The general solution (with Cj fixed to their stationary state values) is a superposition of all 
modes, i.e. 

k 

A(i)=ffJ.(f) = ^/3,e-2-'-*, (83) 

r=l 

^ 1 -2iCrt 

Note that these equations contain only Cr such that Im(cr) < 0, same as in Eq. (j77p . to insure that 
the deviations are indeed small at large negative t. Comparing Eqs. (j83|) and (|77|) . we find 

/3^= ^^"-^"-)n.(cr.-c:) ^,^^,^^^ lm(c^)<0. (85) 

Eqs. (ISSp and ()84p specify deviations of spins from their normal stationary state positions necessary 
to generate the /c-normal-soliton (|75p . Indeed, an arbitrary choice of real Um, 4>m, and large negative 
t = to determines /3m and deviations of spins (f8^ . Equations of motion (f20|) started at t = to with 
these initial conditions produce the /c-soliton solution ()75p with the same values of Cj as those in 
the stationary state. On the other hand, note that generic deviations of spins will modify Cj, see 
e.g. the text following Eq. ()39p . and will not lead to solitons. 

C. Examples of 1 and 2-normal-solitons 

In this subsection, we consider k = 1 and k = 2 normal solitons in more detail, see also the 
Introduction. 

1-normal-soliton. The single normal soliton solution (Fig. [1]) has been previously found in 



Ref. |27|. Here we derive it from the general fc-soliton (|75p as its simplest particular case to illustrate 
our construction of multi-soliton solutions. In this case k = 1 and L^(u) has two complex double 
zeros ci = C2 = /U + ^7 as illustrated in Fig. [9l The corresponding normal stationary state has a 
single discontinuity in the z component of spin (inset in Fig. [9|), i.e. it is the Fermi ground state, 
see Sec. IIIIBI We have seen that in the particle-hole symmetric case 2s| = — sgnej, fi = 0, and 
7 = Ao/2. 

Eqs. (USD, dn]), and dMD yield 

A{ci) = e"+^'^, A{c2) = -6-"+*"^, F = -i —^ 2it,t-i^^ 

cosh(27t + q) 



and 



\A{t)\ = 2\F{t)\ = -^ -. (86) 

' ^ ^' ' ^ ^' cosh(27t + a) ^ ^ 
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Graphically, the single soliton is represented by a single peak located at to = —0/27, see Fig. [TJ 
The parameter 7 controls the width and the height of the peak. 

There is 2A; — 1 = 1 unfrozen separation variable ui. Therefore, Sk{u) = u — ui and rfc_i(u) = 1. 
Eq. ([68]) implies 

ui (t) = fj. — Z7 tanh(27t + a) . (87) 

Note that ui -^ fizizi^ = ci^2 as t — > =Foo in agreement with the results of the previous subsection. 
The separation variable starts from the complex zero ^ + 27 of L^(n) at t = —00 and goes to the 
complex conjugate zero // — Z7 at t = 00 along the straight line connecting the two zeros shown in 
Fig.d 

Individual spin components can be derived from Eqs. (i56]) . (f59|) . and (f87l) . We have 



1 (n-Ci)(n-cJ)i?„_2H r / \ T (^i-'"l)-Rn-2('u) 

Therefore, 

L^{u) = -A{t) ^^"^^^ L„(n). 

(u-ci)(u-c5;) 

Using expression (j4ip with /j = — sgne, and comparing the residues at poles at n = €j on both 
sides of the above equation, we obtain[27] 



.7(t) = 4(t) + i4(t) = A(t) 



[ej-ui{t)]sgn tj 



2[(6,-/.)2 + 72]- 
Similarly, the second equation in Eq. ([58]) yields 

|A(t)P ^■ 



sm = fi^ 



(e, - //)2 + 72 

2-nor'mal-soliton. Now k = 2 and L2(n) has four complex zeros. Fig. I10[ The limiting excited 
normal state exhibits 2A; — 1 = 3 jumps in s^(ej), see the inset in Fig.[Tni In Sec. lIIIBl we considered 
such a stationary state with 2s j = — sgnej(e^ — a^) and determined the corresponding complex 
zeros. 

For a < Ao/4 these zeros are purely imaginary (Fig. [TOb). ci = i'yi, C2 = ^72, C3 = —i'yi, and 
C4 = —^72, where 71^2 are given by Eq. ([55]) . Eq. ([75]) yields 

h{t) 



\A{t)\=A 
where A = A\^2 — 7i| and 



h{t)h{t)-h?{t) 



(88) 



^^^^ ^ gi0, cosh(27it + ai) ^ ^^^^ cosh(272t + 02) _ ^^^^ 

271 272 
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The plot of the 2-normal-sohton (j88p displays two peaks, see Fig. [2l Parameters ai^2 determine 
the location of the peaks in time, while 71^2 control their widths and heights. The 2-soliton can be 
viewed as a nonlinear superposition of two single solitons. At large separation between solitons in 
time, |ai — 02] ^ 1, we obtain from Eq. 



|A(t)| Rs — \ — , (90) 

cosh(27it + ai + r/) cosh(272t + a2 — ??) ' 

where the phase shift r] is 

tanhr/ = sgn(Q2 — aij^j ^. 

7i + 72 

In deriving Eq. (j90p we neglected the terms of relative smallness e"!"^""^!. We see that at large 
separation, the 2-normal-soliton reduces to a simple sum of two single solitons as shown in Fig. [2l 
This is a general property of solitons and one can show that the general /c-normal-soliton ()75p also 
obeys this rule, see e.g. Fig. [3] and Eq. (112p . For small separation the two peaks merge into one. 
When a > Ao/4 in Eq. (j55p the four roots of L2(u) have the form it^ it i'y (Fig. llOb). where 



^ = Ao/4 and 7 = y^a^ — Aq/16. In this case the 2-normal soliton is again given by Eq. 



where now A = IG/iy^/x^ + 7^ and 



h(t) = e-2*/^*+*<^i cos'^(27^ + Qi - il3) ^ ^2iMi+i<^2 cosh(27t + 02 + iP) .g^s 



An additional feature as compared to Eq. (j89p is that here the two terms "rotate" with frequency 
4/i with respect to one another. For large separation, \ai — 02! S> 1, this has no effect - the plot 
of |A(i)| still shows two peaks well separated in time, dashed lines in Fig. [2)o. Now the peaks are 
the same, i.e. 71 = 72 = 7 in Eq. (j89p . In contrast, when the separation is small there is a single 
peak as in the 2-soliton (j89p but with an amplitude modulated by an oscillation with frequency 
a; ~ 4/x = Aq, see Fig [2b. 

V. ANOMALOUS SOLITONS 

In this section, we construct 1- and 2-anomalous-solitons (see also Figs. |3]and[5]) - solutions of 
Bogoliubov-de Gennes equations for |A(t)| that asymptote to anomalous stationary states (j23p 
as t ^ ±00. These solutions show the same solitonic signatures as normal solitons, see the 
introductory paragraph in Sec. IIVI In particular, they are expressed in terms of exponentials 
and multi-solitons break up into a sum of well separated single anomalous solitons in a certain 
limit. 
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A. Single anomalous soliton as a special case of a 3-spin solution 

A single soliton corresponds to the anomalous state with one unstable mode in the linear 
analysis, i.e. L^('u) has two double complex zeros in addition to single zeros u = izAa, see 
Sec. IIII Al and Fig. [71 We considered a state of this type in Sec. IIIIBI In this example, spins in the 
energy interval (—0,0) are flipped; e^ = sgn(|em| — a) in Eq. ([23]) as shown in Fig. [71 (inset). In 
other words. Cooper pairs for single particle states —a < e < a are excited. This state is particle- 
hole symmetric ()44p and the complex zeros of L^('u) are therefore purely imaginary, u = ±17. As 
we have shown in Sec. IIII A| this anomalous state is unstable for 7 > A^. 

For the particle-hole symmetric case equations of motion (j20p have the following form: 

sj = -2e,sj, s,^ = -2AsJ, sj = 2A4 + 2e,sJ, (92) 

where A = g^ ■ sj is real since ^ • s^ = at all times. Let us solve Eq. ([92|) under the condition 
that at t ^ —00 the solution asyrnptotes to the above anomalous state. As mentioned below 



Eq. (|32p and detailed in Refs.l29landl33l. when L^(u) has m complex conjugate zeros (the remaining 
2n — 2m zeros are real) the problem is reduced to solving equations of motion (j92p for m effective 
spins. In the present case m = 3 (counting the pair of double zeros as two pairs) and therefore we 
will need to solve Eq. (|92p for three spins. 

This reduction can be seen in Eqs. ([M]) and ([35]) . Suppose Q2n{u) has only three pairs of 
complex conjugate roots (ci, c^), (c2, Cj), and (03, Cg). There are only 3 — 1 = 2 unfrozen separation 
variables, while the remaining n — 3 are frozen into the n — 3 double real roots of Q2niu), see the 
text following Eq. ([32]) . Suppose c is a real root and let u„_i = c. Then Q2n{uj) contains a factor 
{uj — c)^ which cancels Uj — Un-i = Uj — c in the denominator of Eq. (^^. This cancellation occurs 
for all frozen separation variables and we obtain 

j_ = 2iJ_(ui + n2), (94) 

where Qe{u) = ni=i(^ ~ Ci){u — c*). Eq. ([M]) follows from Eq. ([55]) . since Y^j^j, Jz, and the 
sum of frozen separation variables, X]?=3 % vanish due to the particle-hole symmetry [53]. We see 
that equations of motion (|93p and (|94p are exactly the same as ([34p and (|35p for n = 3 in the 
particle-hole symmetric case. Since the latter equations and Eq. ([92]) are equivalent, Eqs. ([93]) and 
([94]) describe the motion of three effective spins Si, S2, and S3. Note that J-(t) and consequently 



A(i) = gj-{t) are the same in both problems. Moreover, one can show 29l. l33| that the original 
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spins are linearly related to the effective ones, i.e. 

Sj = OjSi + bjS2 + djSs (95) 

Thus, to construct a single anomalous soliton, we need to solve Eq. (|92p for three spins. 

First, let us obtain a general 3-spin solution for which Qq{u) has three distinct pairs of complex 
conjugate roots. As discussed above, the soliton corresponds to the special case when two of these 
pairs, ±^7, are degenerate. The third pair is u = ±iAa and therefore Qq{u) = {v? + 7^)^(n^ + A^). 
The particle-hole symmetry of the 3-spin problem implies ei = — e, e2 = 0, es = e and 

ox _ qx — q qy,Z _ qy,Z _ q qX _ ^ qy,Z _ p. /Qf>\ 

'-'l — '^3 — '^xj '-'i — 1^^ — '-'y,z^ '-'2 — 7T' '-^2 — v^^J 

Using A = gJ2m=i ^rn = '^d^x — 5/2 and integrating Eq. ([92]) . we determine the effective spins 

A 1 A A2 



where C is an integration constant. Combining Eqs. (j97p and (j95p . we derive the original spins in 
terms of A(t), 

s^ = AjA + Fj, sy = BjA, s'j=CjA'^ + Dj, (98) 

where Aj, Bj, Cj, Dj, and Fj are time-independent. The constants Bj, Cj, and Dj are odd in Cj, 
while Aj and Fj are even by particle-hole symmetry ()44p . i.e. Bj = B[ej) = —B{—ej) etc. Since 
A = qY^-s^ we also have 

n n 



Eq. (j98p is similar to the ansatz of Ref. I27l . which is obtained by setting Fj = 0. Nevertheless, this 
difference is important as this ansatz yields 2-spin solutions [29], while here we construct 3-spin 
ones. 

Substituting Eq. ([98]) into equations of motion ([92|) , we find 



A, = -2e,B,, C,=-B„ F, = ^^^, D, = 2{e] - C2)B, 



Bj = - lJ- =, Qe{u) = u\u^ - C2f + 4 - C3U^ 
4vQ6(ei) 



(100) 



where ej = ±1. Since Bj is odd, Sj must be even, e(ej) = e(— ej). The polynomial Qg{u) is 
the same spectral polynomial that appears in Eq. ([93|) . see the discussion of few spin solutions in 
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Refs. |29| and l33l. The coefficients of the polynomial Qq{u) are constrained by Eq. ([99]) . Plugging 



Eq. (fTUnil into Eq. i^, we obtain 



e?eo 



y^=^ = o, y-x^ = -, (101) 

which provides two constraints on three parameters ci, C2, and C3. Thus, 3-spin solutions con- 
structed here are a one parameter family of solutions to Eq. (|92p . 

It remains to determine A(t) for 3-spin solutions. The equation for A(t) can be obtained from 
the condition that the length of spins is conserved by the evolution, s^ = 1/4. With the help of 
Eqs. ([98]) and (|100p this condition reduces to 

A2 = -P4(A), P4(A) = A^ + 4c2A2-8ciA + 4c3. (102) 

For general -P4(A) the solution of this equation is an elliptic function. Here we are only interested 
in an anomalous soliton. As discussed above, it corresponds to a special choice of the spectral 
polynomial Qq{u) = (n^ -|- 7^)^(u^ -|- A^). According to the expression for Qq{u) in Eq. (|100p . this 
implies 

2a a 2 ^a 2a2 

Cl = -7 Aa, C2 = ^-7, C3 = — -7A„. 

For these values of the parameters, Eq. ([102p for the order parameter takes the form 

A^ = -(A - A„)2(A2 + 2AAa + A^ - 47^). (103) 

Now the fourth order polynomial on the right hand side has a double root A^, which means that 
Eq. ([103p can be solved by elementary means. Note also that the stationary state value A(t) = A^ 
is also a solution. In terms of a new variable y = (A — A^)^^ Eq. (|lU3p reads ij^ = A^y^ — 4Aay — 1, 
where A = 1\J"{^ — A^. We obtain 



A 



2 



A(t) - Aa = — — -, A = 2V72 - A2. (104) 

^ ' 2Aa ± 27 cosh(At + a) ' ^ ' " ^ ^ 

The constraints (jlOip become the gap equation ([24p and the equation determining imaginary zeros 
±i7 of L2(u) 

e,- 2 



^— ' / O 0\ /~T I~0 ' ^— ' 



, (e,2 + 7^)^e^2 + A^ TVS^ + ^^ ^ 

We solved these equations in Sec. IIIIB] see Eqs. ([18]) and ([50]) . 
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Finally, Eq. (j98p together with Eqs. (jlOOp and (jl04p yield the individual spin components for a 
single anomalous soliton 



2(eK7^)J.KA2 2^e^2 + A2 



.s//,^ _ gjgj^(0 



'j W = , r^ => (105) 



^.(^) _ e,e,iA^it)-Al) e,e, 



Eqs. p04p and (jlOSp describe a single anomalous soliton solution to the equations of motion ([T|) and 
(j20p . Note that for i — > ±00 the order parameter A(t) — > A^ and the spin components tend to their 
stationary state values (j23p . For A^ = the anomalous soliton (jl04p turns into the normal one 
(|86p . Graphically, the anomalous soliton is represented by a single peak similarly to the normal 
soliton, see Fig. HI Parameters Aq and 7 control its width and height, while a determines its 
position in time. 

B. 2-anomalous-soliton solutions 

Two and higher anomalous solitons can be derived by solving equations of motion for the 
separation variables p6p similarly to the construction of normal solitons above. The fc-anomalous- 
soliton corresponds to a root diagram of L^(n) with k double complex zeros ci, . . . ,Cfc and a pair 



of single zeros HAa- Then, the denominator of Eq. (|36p is y^u^ + A^ Yl^{u — Ci), i.e. only a 
second order polynomial remains under the square root. In this case, Eq. (j36p can be integrated in 
elementary functions. However, here we restrict ourselves to the 2-soliton case and adopt a simpler 
approach to construct it. 

Spin components of the anomalous stationary state, to which the fe-soliton asymptotes at large 
times, display 2k discontinuities. We analyzed an example with four discontinuities in Sec. IIIIBl 
see also Fig. [51 In this example spins in energy intervals (—6, —a) and (a, 6) are flipped, which 
means Cm = sgn(|em| — a)(|em| — b) in Eq. (p3]) . L^(ii) has four complex double zeros ±i^i^2 in 
addition to single zeros iiA^ as illustrated in Fig. [8l We assume 72 > 71 > A^. Therefore, there 
are two unstable modes in linear analysis with growth rates A1.2 = 2(7^2 ~ ^a)^ j ^^^ Sec. IIIlAl 
The 2-anomalous soliton must have the following properties: a) A(t) -^ A^ as t — > ±00 while 
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at large t it should reproduce the linear analysis, b) for Aq = it should be equivalent to the 
2-nornial-soliton described by Eqs. (j88p and (|89|) . and c) in a certain regime the 2-soliton must 
break up into a sum of two single solitons ()104p . This suggests the following ansatz for the 2-soliton 

A{t) - A, = jjlj^, (106) 

f = 0-0 + — cosh(Ait + ai) + — cosh(A2t + 02), (107) 

Ai A2 

where qq, ai, and 02 are time-independent parameters. To determine them, we require that for 
1 02 — Oil ^ 1 the 2-soliton be well approximated by a sum of two single anomalous solitons (cf. 
Eq. dHOD) 

\2 \2 

A(t] — A ~ — \ — flOS") 

^' " ^ 2Aa ± 271 cosh(Ait + ai + r?) 2Aa ± 272 cosh(A2t + 02 - r?) ^ ' 

Neglecting terms of relative smallness e"!"^""^! in Eq. (|106p . we indeed obtain Eq. (jlOSp when 
tanh(ry/2) = A1/A2 and 

2Aa 271 ,^ , 272 

T2T2 ±72772 r2TCOsh(Ait + ai)± 

^1^2 '^ll^2~'^lJ ^2 1^2 ~ ^2 J 



/ = T27I =^ T27T2 r2TCOsh(Ait-Fai)ifc ^ ^ -27- cosh(A2t + a2), Ai,2 = 2^/7^3 " ^a- (109) 



Further, one can verify that the 2-anomalous-soliton given by Eqs. (jl06p and (jl09p also has the 
properties a) and b) discussed above. Its plot consists of two peaks levelling off to the stationary 
value Afl at large times, see Fig. [5j The amplitudes and the widths of these peaks are determined 
by parameters 71, 72, and A^. 

VI. CONCLUSION 

In this paper, we constructed soliton solutions of time-dependent Bogoliubov-de Gennes equa- 
tions ([1]) or, equivalently, Gorkov equations (|20p describing the collisionless dynamics of a fermionic 
superfluid. There are two types of solitons. Normal solitons asymptote at t ^ ±00 to normal sta- 
tionary states (j27p . which are simultaneous eigenstates of the mean-field BCS Hamiltonian ([3|) and 
the Fermi gas. These states are characterized by zero order parameter, A = 0. We have derived 
the general /c-normal-soliton solution, Eqs. (|7Up . (|7ip . and (|75p . and matched the soliton constants 
to small deviations from the corresponding stationary state. We considered the 2-soliton example 
(jSSp in detail and related its parameters to those of the asymptotic stationary state. Examples of 
A; = 1, 2, and 3 normal soliton solutions are shown in Fig. [H [51 and El At large separation between 
the solitons, the /c-soliton becomes a simple sum of k single solitons, see e.g. Eq. (j90p and Figs. [2] 
andO 
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Anomalous solitons asymptote to unstable eigenstates of the mean-field BCS Hamiltonian (j23p 
with nonzero anomalous average. We have obtained one, Eqs. pU4p . (|1U5|) . and Fig. [H and two, 
Eqs. (jl06p . (jl09p . and Fig. [5l anomalous soliton solutions and related their parameters to those 
of the corresponding stationary states. The single soliton is a special case of a more general 3- 
spin solution, which we have also derived. In the vicinity of a stationary state, both normal and 
anomalous multi-solitons break up into a sum of single solitons. These single solitons are unstable 
normal modes in the linear analysis around the stationary state. 

The utility of the soliton solutions is that they are explicit and are in terms of elemen tary 



functions (exponents), in contrast to the general solution in terms of hyperelliptic functions [29|]. 
At the same time, the dynamics in many physical situations is multi-soliton. The combination 
of these two factors makes solitons potentially quite useful in various problems in non-stationary 
superfluidity. Consider, for example, the collisionless dynamics triggered by an abrupt change 
of the pairing strength. In most cases of interest L^(m) has only few isolated zeros, while the 
remaining complex zeros merge into continuous lines [3^]. We believe that the latter zeros can 



be treated as being degenerate and their contribution is therefore mu 
then a superposition of a (quasi-)periodic few spin solution 



291, 



30, 



ti-soliton. The solution is 



33l | with a multi-soliton one. 



Superpositions of this type are referred to as solitons on a (quasi-)periodic background in the soliton 



theory 



53( 1 . In particular, when the system is in the ground state before the coupling change, the 



collisionless dynamics governed by Eg. ( ^ can produce asymptotic states with a constant nonzero 



order parameter or a gapless state 



M, 



35|, 



36| |. In these cases L^(u) has either a single pair of 
nondegenerate zeros or no such zeros. Therefore, according to the above reasoning, the dynamics 
leading to these asymptotic states is described by a multi-soliton solution of a normal type for the 
gapless state and of an anomalous type otherwise. 
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